Large-Scale Structure in Brane-Induced Gravity 
II. Numerical Simulations 
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We use A''-body simulations to study the nonlinear structure formation in brane-induced gravity, 
developing a new method that requires alternate use of Fast Fourier Transforms and relaxation. 
This enables us to compute the nonlinear matter power spectrum and bispectrum, the halo mass 
function, and the halo bias. From the simulation results, we confirm the expectations based on 
analytic arguments that the Vainshtein mechanism does operate as anticipated, with the density 
power spectrum approaching that of standard gravity within a modified background evolution in 
the nonlinear regime. The transition is very broad and there is no well defined Vainshtein scale, 
but roughly this corresponds to fc« ~ 2/iMpc~^ at redshift z — 1 and fc, ~ IfeMpc"'^ at z = 0. 
We checked that while extrinsic curvature fluctuations go nonlinear, and the dynamics of the brane- 
bending mode C receives important nonlinear corrections, this mode does get suppressed compared 
to density perturbations, effectively decoupling from the standard gravity sector. At the same time, 
there is no violation of the weak field limit for metric perturbations associated with C. We find 
good agreement between our measurements and the predictions for the nonlinear power spectrum 
presented in paper f, that rely on a renormalization of the linear spectrum due to nonlinearities in 
the modified gravity sector. A similar prediction for the mass function shows the right trends. Our 
simulations also confirm the induced change in the bispectrum configuration dependence predicted 
in paper 1. 



I. INTRODUCTION 

There is strong observational evidence from different 
techniques that the universe is currently undergoing ac- 
celeration. The most studied explanation of cosmic ac- 
celeration is that the universe is filled with an addi- 
tional stress-energy component, usually called dark en- 
ergy. The simplest example is the cosmological constant, 
while more general possibilities include scalar fields with 
suitable potentials, as in early universe inflation. Since 
the study of the cosmology relies crucially on Einstein's 
general relativity (GR), and we test GR stringently only 
up to the scale of the solar system, it is possible that 
deviations from GR occur at cosmological scales. The 
observed cosmic acceleration may well be due to the fact 
that GR does not hold at such scales. It is thus important 
to explore predictions for theories that modify gravity at 
large scales to explain cosmic acceleration. 

Among these modified gravity theories, brane- 
induced gravity in five dimensions, the so-called Dvali- 
Gabadadze-Porrati (DGP) model [1] is one of the most 
widely studied. In this model, gravity is intrinsically 
five-dimensional, and a standard four-dimensional Ein- 
stein term is induced by the presence of standard model 
particles in a four-dimensional brane. Thus, while at 
small scales gravity reduces to GR, at large scales gravi- 
tons can leak out to the extra-dimension and the gravita- 
tional force law becomes five-dimensional. Interestingly, 
although it was not put forward to explain acceleration, it 
was later found that this model exhibits self-acceleration 
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in the late universe [2] . This makes this model very inter- 
esting phenomenologically as it may explain the observed 
cosmic acceleration without invoking dark energy. In this 
paper we focus on the DGP model, but it is expected 
that the techniques developed here will be applicable to 
other models of massive gravity. For a review of the DGP 
model, see e.g. [3-5]. 

Distinguishing dark energy from modified gravity is 
thus one of the main questions in understanding the cause 
of cosmic acceleration. At the level of expansion history 
alone, it is very difficult to do so, as any modified gravity 
expansion history can be mimicked by a suitable equa- 
tion of state for the dark energy. However, the growth 
of perturbations will typically be different [6] , unless one 
tunes the properties of the dark energy to make it cluster 
(instead of being smooth) on subhorizon scales to mimic 
a modified gravity model (see e.g. [7] for such an attempt 
in linear perturbation theory); however, such fine tuning 
may be even more contrived at the nonlinear level. In 
a companion paper ([8], hereafter paper I), we discuss 
predictions at the nonlinear level based on perturbation 
theory (see also [9] for closely related work), here we ex- 
tend the work in paper I by running cosmological TV-body 
simulations. 

Studies of structure formation using A^-body simula- 
tions in modified gravity has been carried out in [10-12], 
but were restricted to linearizing the Poisson equation, 
while the field equations of viable modified gravity theo- 
ries are typically nonlinear, since it is through nonlineari- 
ties how the extra degrees of freedom that modify gravity 
at large scales get suppressed at small scales, recovering 
GR as required by local tests. In [13-15], A^-body simu- 
lations are carried out for the /(i?) models, which involve 
nonlinear couplings of the scalar field to matter analogous 
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(though very different in detail) to the second-derivative 
self-interactions present in the DGP model. 

In this work, we concentrate on fully nonlinear TV-body 
simulations of the DGP model, solving numerically the 
equations of motion derived in paper I [8] (see also [16]). 
The growth of perturbations in the linear and nonlinear 
case was first studied using spherical symmetry in [6]. 
The linear perturbations were studied beyond spherical 
symmetry by [17], where the behavior in the bulk was 
established. All these treatments use the quasistatic ap- 
proximation, appropriate for scales below the Hubble ra- 
dius. For larger wavelengths, numerical solutions of lin- 
ear theory were derived in [18, 19]. 

One of the interests of the present work is to study 
in some detail the so-called Vainshtein mechanism [21] 
in the cosmological context. In paper I we studied it 
analytically, here we resort to A'^-body simulations for a 
solution of the fully nonlinear equations that describe it. 
The Vainshtein mechanism was originally found in mas- 
sive gravity. In the linearized case, the graviton prop- 
agator in the limit that the graviton mass goes to zero 
differs from the massless graviton theory: this is the well 
known vDVZ discontinuity [20]. As conjectured by Vain- 
shtein [21], however, the discontinuity should disappear 
in the exact non-perturbative solution when nonlincar- 
ities are included. While there have been some doubts 
in literature about whether the Vainshtein mechanism 
works in massive gravity (see e.g. [22]), recent work con- 
firms its validity and clarifies the subtleties involved in 
matching the expected small-scale and large-scale behav- 
iors [23]. 

The interest in the framework of cosmic acceleration 
is that this mechanism is responsible for modified grav- 
ity theories becoming GR at small scales, and the tran- 
sition scale is important to pin down quantitatively to 
know where we can best test modifications of gravity in 
observations. In the DGP model, the Vainshtein effect 
is naturally built in [24]. For a compact object, within 
its Vainshtein radius = (r^Tgch)^^'^, where rgch is the 
Schwarzschild radius and Tc is the crossover scale, the 
extra scalar degree of freedom is strongly coupled and 
heavily suppressed, so that massless gravity is recovered 
in the nonlinear regime [24-28] . In the case of cosmology, 
the precise value of is not as easy to obtain since it 
depends on the details of the fluctuation spectrum, while 
analytic estimates exist [6, 8, 16] its accurate determina- 
tion requires a numerical simulation to cross check the 
analytic estimates. 

Our aim here is not to make precise predictions for the 
DGP model, as it is challenged from observations [30-32] 
(see however [33] for recent results with different conclu- 
sions) and theoretical considerations [29], but rather un- 
derstand the nonlinearities in the modified gravity sector 
in this model, since there are very good reasons to expect 
that it is a representative example of more complicated 
theories of massive gravity that may have better theoret- 
ical and observational properties, e.g. [34-36]. In partic- 
ular, it is expected that the method introduced here for 



solving the nonlinear field equations would be useful in 
such models. 

The rest of the paper is organized as follows. We first 
recall the background and perturbation equations for the 
nonlinear DGP model in Sec. II. In Sec. Ill A, we out- 
line the general procedures for A^-body simulations in 
modified gravity. The FFT-relaxation method is first de- 
scribed in Sec. IIIB and then generalized in Sec. IIIC. 
A brief description of the computer resources require- 
ment for the FFT-relaxation is also given at Sec. HID. 
We present the power spectrum, bispectrum, halo mass 
function and halo bias from our simulations in Sec. IV, V, 
VI A, and VI B respectively. We conclude in Sec. VII. For 
more discussion on technical aspects, in particular the el- 
lipticity of the field equation, and the comparison of the 
FFT-relaxation method with the standard Gauss-Seidel- 
relaxation (GS-rclaxation) method, see Appendix A. 



II. EQUATIONS OF MOTION 

The background Fricdmann equation in the DGP 
model in the self-accelerated branch is [2] 
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where pm is the energy density of matter, and is a 
cross-over scale which characterizes the length scale at 
which gravity goes from 4D to 5D behavior. In the long 
time limit, when matter density is sufficiently diluted, 
the universe enters a de Sitter phase with H = r~^. To 
explain the cosmic acceleration, has to be of the order 
of the present Hubble radius. In particular, one can write 
the current Hubble rate Hq in terms of the present matter 
density parameter fl'^ as 
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The evolution of scalar metric perturbations in the sub- 
horizon limit can be reduced to a coupled system of equa- 
tions involving the Newtonian potential (f> and a metric 
potential C that characterizes the extrinsic curvature of 
the brane in the subhorizon limit. The modified Poisson 
equation reads [8] 
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where 5 is the density contrast and C obeys a fully non- 
linear partial differential equation 



2r/ - 1 



3(77-1) 



(l-/3V-i)5, 



(4) 



3 



see paper I for a full derivation. Here wc have defined 

1] = rcH, (5) 
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and the dimcnsionless operator 
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has large eigenvalues {k/aH ^ 1) in the subhorizon limit. 
The nonlocal operators V — and = 1/V— can 
be easily handled as they are simply k and k~^ in Fourier 
space. The major obstacle in solving Eq. (4) stems from 
the nonlinear terms (V^C)^ and (VijC)^. One of the 
main aims of this paper is to present a convergent algo- 
rithm to solve Eqs. (3,4). 

It will also be useful to linearize the field equations 
since we shall compare the nonlinear solution against the 
result of linearizing the field equations, but keeping the 
equations of motion for matter fully nonlinear. After 
dropping the nonlinear terms, we can express V^C in 
terms of the density contrast 5. We can then plug in 
the expression for V^C into Eq. (3), and expanding the 
denominator in series of k~^ up to first order, we get the 
linearized Poisson equation 
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where GcS is an effective scale and time dependent New- 
ton's constant 
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and q is the acceleration parameter 
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We emphasize that we can rigorously absorb the effects 
of modified gravity into a scale and time dependent effec- 
tive gravitational constant only if we drop the nonlinear 
terms. When nonlinear effects are included, the argu- 
ments of paper I suggest that one can represent most 
(though not all) of the nonlinear effects by a renormal- 
ization of G, part of the purpose of this paper is to verify 
those predictions against fully nonlinear numerical sim- 
ulations. The first term in Eq. (9) represents the change 
in G due to the linearized dynamics of C (which medi- 
ates a repulsive force), while the second term inducing 
scale dependence represents the change in the gravita- 
tional force-law as one approaches the cross-over scale 
Tc- Both of these effects make G^s less than G. When 
the Vainshtein effect sets in due to nonlinearities in the 
dynamics of C, Gcs/G is driven to unity and one recovers 
the standard Poisson equation 



but the expansion history will still be modified. There- 
fore, the benchmark to study the Vainshtein effect is to 
check whether the growth of perturbations at small scales 
approaches that of standard gravity with a modified ex- 
pansion history. 



III. iV-BODY SIMULATIONS 

We shall first briefiy review the general method of par- 
ticle mesh (PM) iV-body simulations in Sec. Ill A. We 
then introduce and describe the FFT-relaxation method 
we use in this paper in sections III B-III D. For a compar- 
ison between this method and the standard Gauss-Seidel 
relaxation method see Appendix A. 



A. PM iV-body Simulations 

Since the procedure for running modified gravity iV- 
body simulations is largely similar to the standard GR 
one, which is well-known, see e.g. [37-39], we will only 
briefly outline it here. We also state the value of the 
cosmological and internal parameters used in our code. 

Our simulations use 256^ particles and start at z = 49 
for 128 Mpc/h box size and z = 24 for the larger 
boxes. The initial conditions were generated using sec- 
ond order Lagrangian perturbation theory [40] with the 
transfer function from CMBFAST [41]. Unless otherwise 
noted our simulations assume f2„i = 0.27, h = 0.7 and 
rif, = 0.04. The normalization is set by the present lin- 
ear rms fluctuation erg, and we use erg = 0.9 for the case 
of standard gravity, with the initial amplitudes the same 
for the different runs (thus erg = 0.715 at z = for DGP 
models). 

The computational domain is a cube with periodic 
boundary conditions. Due to the limitation of computer 
resources, we typically use 512'^ grid points in our flnal 
results unless otherwise specified. Thus the number of 
particles are 1/8 of the number of grid points. Each time 
step is uniform in In a, where a is the scale factor, of 
width 0.005. When starting from z = 24, it takes 618 
time steps to reach z = 0. 

The method we use is the standard Particle-Mesh 
(PM) algorithm. In each step we interpolate the mass 
to the grid points from the particles in order to calculate 
the potential and hence the force at the grid point. The 
force is then interpolated back to the particle from the 
grid points using the same interpolation scheme. The 
assignment function we use is the Cloud-In-Cell (CIC), 
which is widely used in PM simulations. In ID, the win- 
dow function is 
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M if 1^1 < d 
otherwise 
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where d is the size of the grid. For 3D, the assignment 
function is W {x)W {y)W {z) . Thus the mass distribution 
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at the grid points Pgrid can be given in terms of the par- 
ticle distribution Pparticic as 

Pgiid{x,y, z) ^ J dx'dy'dz' X 

Wix - x')Wiy - y')W{z - z')Pparticic(x', y', z'). (13) 

From the density contrast at the grid points, one can 
calculate the potential by solving the field equation. In 
GR; the field equation is just the Poisson equation, which 
can be solved by a single Fast Fourier Transform (FFT). 
In the DGP model, the field equations are nonlinear, and 
we will develop an algorithm to solve them in the next 
section. 

From the Newtonian potential (j), we can advance the 
particles using Newton's law 



— = -2i7u--V0, (15) 
at 

where x is the comoving distance. The second equa- 
tion holds in modified gravity because in a metric theory 
the motion of non-relativistic particles is sensitive to the 
Newtonian potential as in GR through the geodesic 
equation (with the background spacetime being FRW due 
to homogeneity and isotropy). 

One can calculate V0 from (?!)(x) in real space using 
finite difference implementations of the gradient. But 
we choose instead to compute V(f) in Fourier space. In 
Fourier space, V0 is just zk(/), and we Fourier transform 
ik(/) back to real space to get the force at the grid points. 
Again, we have to interpolate the force at the mesh points 
to the particles. To avoid self-force, we use the same 
assignment scheme as before, thus the force interpolation 
is the same as Eq. (13), except that the roles of grid points 
and particles are interchanged. Once the force has been 
evaluated, we evolve the particle positions and velocities 
using the second-order leapfrog method. 

As a test, in Fig. 1, we compare our code against Gad- 
get2 [42] for the case of GR, where we show the ratio 
of the measured nonlinear power spectrum to the initial 
power spectrum scaled to the corresponding redshift by 
the linear growth factor D^{z) as a function of wavenum- 
ber k. In this test our code is run in a box of size 
1280 /i^^ Mpc a side using 256'^ particles. The power 
spectrum for Gadget has been averaged over 50 realiza- 
tions using different initialization seeds (using 640'^ par- 
ticles in a 1280/i~^Mpc box, see [43] for more details), 
and the associated error bars are one standard deviation 
on the mean of these realizations. The large-scale modu- 
lation of the Gadget result is mainly due to the fact that 
we divided by the average initial power spectrum used 
in our code, which corresponds to only four realizations. 
The range in k shown in Fig. 1 is chosen to go up to half 
of the particle Nyquist frequency, within which the dif- 
ference from Gadget is a systematically lower power by 
less than about 10%. Beyond the range shown in the fig- 
ure, we find that our code gives significantly lower power 



1.3 




10"^ 10"^ 



k (/i/Mpc) 

FIG. 1; Comparison of the measured power spectrum from 
our code against the Gadget2 runs with the same initial con- 
ditions, but of different initial random seeds, for standard 
gravity. The upper panel is for z = 1, and the lower panel 
corresponds to 2; = 0. The Gadget power spectrum has been 
averaged over 49 realizations, while the spectrum from our 
code is from four realizations. The numerical power spectra 
have been divided by the average initial spectrum used by our 
code scaled to the corresponding redshift using linear theory. 

than that from Gadget as expected from smoothing due 
to CIC interpolation. By taking the ratio of the power 
spectrum in Gadget to our code, we see that the ratio 
is constant with redshift, as expected from such smooth- 
ing. This systematic lower power have also been noted 
in many other PM simulations, e.g. [10, 12, 14] when 
making comparisons with high resolution codes or stan- 
dard fitting formulae. In our case, these systematics are 
not too worrisome as we arc interested in understanding 
differences (i.e. ratios) of power between GR, GR with 
modified expansion history (hereafter GRH), and modi- 
fied gravity with linearized and full Poisson equation sim- 
ulations, and for all of them we will use our code. 



B. FFT-Relaxation Method 

In this section, we describe the procedure that we use 
to solve Eqs. (3-4). We shaU first solve Eq. (4) to get V^C 
as a function of 5, and then substitute it into Eq. (3) to 
compute 4). After obtaining C, Eq. (3) can be solved 
using standard FFT techniques. 

To solve Eq. (4), the key idea is to notice that it can 
be regarded as a quadratic equation in V^C, which can 
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be solved for using the quadratic formula in terms of 
(Vy C)^ and other remaining terms treated as sources. 

These sources, e.g. Vy C and v — V^C, depend in a non 
local way on V^C. But these can be evaluated in Fourier 
space and an iterative relaxation method can be estab- 
lished among them. 

Let us now discuss the details of the algorithm. We 
first start from the linearized version of Eq. (4) 



3/3(7? 
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In Fourier space, this equation becomes 
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For notational convenience, we have used the same sym- 
bol for the quantity in real and Fourier space. From this 
Fourier representation, one can obtain C(k) algebraically 
and this is used as our initial trial solution for the first 
time step in the FFT-rclaxation algorithm. Using that 



VyC(k) = -hkjCik), 



(18) 



and Fourier transforming V,jC(k) back to real space 
yields VyC(x). From VyC(x) we can then calculate 
(VijC(x))^, and then we can update V^C(x) by treat- 
ing Eq. (4) as a quadratic equation in V^C(x) using the 
quadratic formula 



with 
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From Eq. (19), we see that the solution to the equation 
is not unique. However, the solution should reduce to 
that of the linearized equation in the linear limit and 
this holds only if we choose the plus sign. It is also not 
a priori clear that the discriminant + AB is always 
positive. We will come back to this issue in Sec. IIIC. 

With the new V^C, we can repeat the procedure to 
find {S/ijCY and the nonlocal term, and therefore B. 
Using the new i?, we can compute V^C again. This 
loop continues until some tolerance is met. That is, our 
relaxation algorithm can be written schematically as 



-a-t- -^4B("-i) 



(21) 



at step n in the iteration. Since each step involves calcu- 
lating FFT's, we call our algorithm FFT-relaxation. 



10 



FIG. 2: Contributions to the discriminant A that do not de- 
pend on a trial solution for C, in voids [5 — —1) as a function 
of z for two choices of splitting the nonlinear terms. The 
top curve {w = 1/3) corresponds to a splitting based on a 
multipolar expansion of the nonlinear kernel, while the bot- 
tom curve {w — 0) denotes the naive splitting presented in 
Sec. IIIB. Note that in the latter case these contributions are 
negative already for z > 1, and thus a good trial solution for 
C is required to avoid A < 0. 



At each step we calculate the residual i?(x), which is 
the difference between the left hand side and right hand 
side of Eq. (4) . A reference standard for comparing the 
size of the error is obtained by building the following 
quantity 



, 3(7^-1) 
77 



\{1-I3^-')S\), 



(22) 



where the angular bracket denotes average over all the 
grid points. We then normalize the average residual {R) 
to e, and monitor the normalized residual {R)/e as the 
iteration progresses. We stop the iteration when the nor- 
malized residual is less than 1%. After solving for V^C, 
we input it into Eq. (3) and solve for (/)(k) in the Fourier 
space and then Fourier transform back ik(/)(k) to get the 
force. 

Except for the first time step in which we use the lin- 
ear solution as the trial solution in the relaxation process, 
after the initial time step we use the previous time step 
nonlinear solution as the seed for the relaxation proce- 
dure. Using the previous nonlinear solution, it takes on 
average about three iterations for the normalized residual 
to decrease to 1%. 



C. Splitting the Nonlinear Terms 

There is no guarantee that the discriminant A = -I- 
AB in Eq. (19) is always positive during time evolution, 
particularly in voids. To see more clearly the source of 
potential problems, let us rewrite Eq. (4) without the 
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nonlocal terms, which we'll assume for simplicity here 
that do not alter things qualitatively (this is confirmed 
by numerical testing) 



3(r;~l) 



5. 



(23) 



It is easy to see that the sum of the two quadratic terms 
is positive definite (see Eq. 29 below), therefore while in 
linear theory the left hand side of Eq. (23) tracks the sign 
of 6 (recall 77 > 1) in the nonlinear case as voids empty 
and (5^—1 with V^C becoming of order unity, the non- 
linear term may take over and this equation may not have 
real solutions, depending on a{r]). That would mean the 
theory is unphysical, or unstable in voids. However, even 
if this is not the case for the true solution, a particular 
implementation of the relaxation algorithm can be sub- 
ject to these type of problems. 

To see this, let as rewrite the discriminant A for 
Eq. (23) 



A. 



(V.,C)^ 



V 
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In Figure 2, we plot the contributions to A that are in- 
dependent of C, i.e. + 12(77 ~ l)r]~^S (which corre- 
sponds to w — 0, with the parameter w to be intro- 
duced shortly), with J = —1 as a fimction of redshift z 
for = 0.2 (the precise value of il^^ does not signifi- 
cantly change this). We see that for z > 1, the value of 
+ 12(77 ~ l)ri~^S is less than for 6 close to -1. This 
can potentially cause problems if the initial guess is not 
good enough and (VyC)^ is too small, making A neg- 
ative. We shall see later that this heuristic argument is 
indeed quite accurate when we discuss Fig. 4. 



To make the sign of A less sensitive to our choice of the 
trial solution we split the nonlinear terms by introducing 
a parameter w that separates how much of the local non- 
linear piece (V^C)^ is solved for immediately and how 
much is left for the iteration process. Specifically, let us 
rewrite Eq. (4) as 



3/3(^-1) 
27^- 1 



3(7,-1) 



= 0. 
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With this simple splitting, our iterative procedure de- 
scribed in Sec. IIIB remains the same except Eqs. (19-20) 
now become 



2(1 -w) 



(26) 



and 



B' 



= (V,jC)2-7«(V2C)2- 

+ fci)(l_;5V-i)5. 
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Neglecting the nonlocal term, the new discriminant re- 
duces to 



4(1 



w) 



(28) 

For w in between and 1, the factor {1 ~ w) reduces 
the magnitude of the term proportional to 6, with the 
desirable effect that the contributions to A that do not 
depend on C are now positive even for completely empty 
voids at all redshifts. Indeed in Figure 2, we also show 
-I- 8(7; — l)ri~^S (corresponding to w = 1/3) with 
6 ~ —I, and we see that this quantity is now positive 
for all z. Of course, now we have introduced a negative 
contribution into A given by — 7«(V^C)^ (which must be. 



since under correct choice of C the result should be in- 
variant under w) but now one expects the contribution 
coming from C terms to be smaller in amplitude because 
of cancellations between the two terms that depend on C 
and thus less sensitive to incorrect guesses coming from 
the trial solution. 

How do we choose the arbitrary parameter wl The 
most physical choice we can think of is w = 1/3, as it is 
motivated by the spherical collapse model. To see this we 
note that in Fourier space, the nonlinear terms in Eq. (4) 
can be written as 



(V2C)2 - (V,,C)2 (k) 



(29) 



<5D(k - ki - k2) kiki [1 - C(ki)C(k2) , 



where ^ = ki • k2 . We can decompose the angular depen- 
dence in multipoles and write I ~ — (2/3) [1 — 7^2 (m)]; 
where 7^2 (m) = (3/i^ — l)/2 is the second Lcgendre poly- 
nomial. Then, 



(V2C)^-(V,,C)2 = ^(V2C)2-^^2 



(30) 



where we have introduced the operator 1^2, such that 
d^k 



V2[A,B] = 



-zkr 



<5D(k-ki -k2) X (31) 



(27r)3 

P2(m) A{\ii)B{\^2)d'kid?k2 
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and in our case 



V2 



\{^^,Cf~\{^^C)\ (32) 



Thus choosing w = 1/3 means that we separate the non- 
Unear terms in a geometric way, and at each iteration 
we solve for the sphericahy averaged field equation treat- 
ing the quadrupole as a source of perturbation. In other 
words, for w = 1/i the C terms in Eq. (28) describe 
the quadrupole of the field C which should be a small 
quantity (compared to the monopole) and thus making 
A more stable against choices of the trial solution for 
C . This decomposition is reminiscent to the approxima- 
tions made in paper I to initialize the resummation of the 
two-point propagator of the modified Poisson equation. 

Another potentially interesting choice of w is that 
which makes the discriminant A = when (5 = — 1, while 
putting the C terms to zero in Eq. (28). In this case w 
is time dependent and given by 



w ~ 1 — a'^ 



12(7,-1)- 



(33) 



Figure 3 shows the normalized residual for different val- 
ues of w. For the purpose of this illustration we have 
taken z = 24 and a box size of 64 Mpc using a 256^ 
grid, although similar results hold for other choices (as 
long as z > 1). In all the the cases shown in Fig. 3, ev- 
erything is kept the same except the value of w. In par- 
ticular, we start from the same linearized solution, and 
therefore the initial normalized residuals are the same 
independent of w. 

Let us denote those grid points at which the discrim- 
inant is negative as bad points, and define the bad frac- 
tion as the ratio of the number of bad points to the total 
number of grid points. For w = 0, in addition to the nor- 
malized residual, we also show the bad fraction in dashed 
lines (the bad fraction is zero for all the other values of 
w shown). When bad points appear, we set the discrim- 
inant in Eq. (26) artificially to zero so we can continue 
the iteration. For w = at the beginning the normalized 
residual (solid line) decreases as we iterate, but it stops 
decreasing at some point, and increases gradually in a 
slightly oscillatory manner. The bad fraction (dashed 
line) also increases gradually but in a strongly oscilla- 
tory fashion. We found that the trend that the normal- 
ized residual decreases at first and then increases slowly 
is typical when the bad fraction is nonzero. 

We also plot w = 0.25, which corresponds to using 
Eq. (33) for the particular ij for z ~ 24, and it converges 
slightly slower than w ~ 0.5 and the physically motivated 
w — 1/3. When w is close to 1, such a.s w = 0.8 shown 
here, the convergence is much slower, although the bad 
fraction is still zero. Although w = 0.5 converges slightly 
faster than w = 1/3, the difference is negligible in the 
range of the normalized residuals we are interested in, 
and w = 1/3 is better motivated theoretically, thus we 
shall use ui = 1/3 in the subsequent computations. 

The saturation of the normalized residual at about 
1.5 X 10~^ is most likely due to the discretization of 
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FIG. 3: Normalized residual of solving the nonlinear field 
equation for u; = 0, 0.25, 1/3, 0.5, 0.8 as a function of the 
number of iterations. For w = 0, we also show the variation 
of the bad fraction (dashed lines), while the bad fraction is 
always zero for all the other values of w. When bad fraction 
is non-zero {w — 0), the normalized residual decreases for the 
first few iterations, and then it bounces back and increases 
gradually. The normalized residuals decreases most rapidly 
for w = 0.5 and 1/3. The saturation at about 10~* is prob- 
ably due to the discretization of the equation, see text for 
discussion. 
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FIG. 4: The bad fraction and normalized residual for w = 
and 1/3 as a function of scale factor a. For w = 0, the 
normalized residual is high when the bad fraction is non-zero 
in the range up to a ~ 0.6, as expected from the arguments 
made in Sec. IIIB. The normalized residual can be easily 
reduced to below 1% when the bad fraction is zero for the 
splitting choice of to = 1/3. 
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Eq. (4). As the grid becomes more refined, the residual 
calculated using a consistent numerical method should 
approach zero. However, this assumes that the fluctua- 
tion at the grid scale remains unchanged as the size of 
the grid increases. We find that the normalized residu- 
als saturate to roughly the same value using 128'^, 256^ 
and 512'^ grids. However, as the size of the grid increases 
the fluctuations at the grid scale increase, since we probe 
smaller scales. When we increase the grid size and artif- 
ically decrease erg to maintain the level of fluctuations at 
the grid scale constant, the saturated normalized residual 
value indeed decreases. We emphasize that in any case 
the saturated value seen in Fig. 3 is much lower than we 
require. 

In Figure 4 we plot the bad fraction and normalized 
residual for w = (top panel) and w = 1/3 (bottom 
panel) during a full simulation from a 0.04 to a = 1 
using 256'^ grid points in a 64/i~^Mpc box. Whenever 
the bad fraction is nonzero we terminate the iteration 
of the relaxation method when the minimum residual is 
reached. For w = 0, the bad fraction is non-zero from 
the beginning a = 0.04 to about a = 0.6, and the min- 
imum normalized residual that can be attained is also 
rather high. Note that the range where the bad fraction 
is nonzero agrees with the simple arguments leading to 
Eq. (24) and displayed by the w = curve in Fig. 2. 
When we set w = 1/3, on the other hand, the bad frac- 
tion is always zero and the normalized residual can be 
maintained below 1% during the whole simulation with 
only a few iterations of the relaxation procedure. 



D. Computational Resources 

The computing time and memory required to perform 
this type of numerical simulations are much higher than 
the standard gravity A'^-body simulations. For example, 
because of non-local terms in Eq. (27), we Fourier trans- 
form S, V^C, (V^C)^, and (V„C)2 to evaluate B' in 
Fourier space. To compute (Vy C(x))^ from V^C(k) we 
have to loop over all k modes and inverse FFT to real 
space six times. Including calculating residuals, we need 
to perform 15 FFTs in one iteration. For each time step 
in the simulation, we need on average three cycles to at- 
tain the desired accuracy. Thus we have to do about 45 
FFTs to solve the field equations for C in one time step, 
while in standard GR one FFT is required. Computing 
the force requires 3 more FFTs in either case. Therefore 
there is a factor of about 48/4 = 12 in terms of cost from 
FFTs alone. 

We have parallelized our code using OpenMP to speed 
up the computations. When the redshift is high, the non- 
linearity is small, therefore we use the linearized DGP 
Poisson equation and switch to fully nonlinear calcula- 
tions at some epoch depending on the size of the box, 
e.g. at a = 0.1 for the 512 h'^ Mpc box. Using 4 CPUs, a 
simulation with 256'^ particles in standard gravity takes 
about 10 hours, while for the nonlinear DGP model it 



takes about 4 days. When using a 512"^ grid, the FFT- 
relaxation code requires about 16 GB of memory. Al- 
though our implementation is not very economical as far 
as computing resources is concerned, and can certainly 
be improved, it is more than enough for our purposes. 

IV. THE POWER SPECTRUM 

We now present results for the power spectrum of the 
density perturbations measured from simulations. We 
would like to compare the power spectrum for three dif- 
ferent models obtained by simulating different perturba- 
tion equations but with the same background evolution 
given by the Friedmann equation Eq. (1), which we de- 
note by 

P„1DGP(A:), PmGp{k), PGRH(fc), (34) 

corresponding to the fully nonlinear DGP model (Eqs. 3- 
4), the hnearizcd DGP model (Eq. 8), and GR (Eq. 11) 
with the same H{z), respectively. Since all three models 
share the same expansion history, the differences in their 
power spectra arise only from the different gravitational 
forces for perturbations. 

At large scales, we expect that PniDGP — -Pidgp as the 
linearized modified Poisson equation, Eq. (8), becomes a 
good approximation. On the other hand, at small scales 
we expect that PniDGP — ^grh as the Vainshtein mecha- 
nism suppresses C relative to (j) and the standard Poisson 
equation becomes a good approximation. Therefore, the 
extra nonlinearities in the modified gravity sector mean 
that there are two models with the same linear spectrum 
but different nonlinear spectrum (PniDGP and Pidgp) and 
two models with different linear spectrum that have the 
same nonlinear spectrum (PniDGP and Pgrh)- 

Figures 5 and 6 show the power spectra for these three 
models at redshift z ~ 1,0, respectively, with the left 
panels showing power spectra, and the right panels show- 
ing ratios to better see the details. Each set of sym- 
bols corresponds to merging two computational boxes, of 
size 512 /i"^ Mpc (large scales up to A: ~ 0.3 /i Mpc~^) 
and 128/1-1 Mpc (fc ~ 0.3hMpc~^ and up). While 
one should not take the absolute value of the power 
spectra too seriously at high frequencies (since suppres- 
sion due to grid smoothing is significant) we can trust 
much more the relative differences between the spec- 
tra. To see this better, the right panels show the ra- 
tios PniDGP /Pidgp and PcRn/PniDGP- We see the ex- 
pected behavior, PniDGP /Pidgp becomes unity at large 
scales, and gets amplified at small scales as the Vain- 
shtein effect makes gravity stronger. On the other hand, 
PGRn/PniDGP is larger than unity at large scales, as DGP 
gravity is weaker than GR, and the ratio is driven to unity 
at small scales once the Vainshtein mechanism sets in. 
We see that the transition to GR is not complete at the 
smallest scales shown, even at z = and k ~ 6hMpc~^ 
there is still about 10% enhacement. This is in part due 
to the fact that at z = the Pgrh /PniDGP ratio is larger 
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FIG. 5: Dark matter power spectra from the nonlinear DGP model (nlDGP) , linear DGP (IDGP), and GR perturbations with 
the same expansion history (GRH) at z — 1. The left panels show the power spectra, and the right panels shows ratios to 
better see the differences. Two sets of computational boxes are shown for each case, covering a different range in k (see text). 
The solid line denotes the predictions from paper I for PniDGP (left panel) and -PcRH/f'niDGP (right panel). 




k [h/Mpc] k [h/Mpc] 

FIG. 6; Same as Fig. 5 but for z; = 



than at z = 1 at large scales, so the Vainshtein effect may be defined by the scale at which PgrhZ-PhIDGP starts 
has to overcome a larger difference at z = 0. A Vain- to decrease, this is about fc* ~ 2ft,Mpc~^ at z = 1 and 
shtein scale (analogous to in the Schwarzshild case) /c* ~ l/iMpc~^ at z = 0. Note that at intermediate 
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scales the PgrhZ-PhIDGP goes up slightly, this is expected 
as the Pgrh is enhanced by nonlinearities more than 
-PniDGP due to its larger growth factor or erg. 

The solid lines in Figs. 5 and 6 denote the predictions 
from paper I for PniDGP (left panel) and PgrhZ-PhIDGP 
(right panel). We see very good agreement, at the few 
percent level for the ratios, and slightly worse for the 
absolute power (about 5,3% for z = 0, 1, respectively, at 
quasilinear scales), until smoothing due to the CIC kernel 
substantially cuts the power at high-fc. The prediction 
tends to underestimate the N-body simulations, and this 
is expected to some extent from the details described in 
paper I. First, we neglected the extra mode-coupling in- 
duced in modified gravity, which was estimated in paper I 
to be about 1-2% from the bispectrum results. Second, 
the predictions use the GR HaloFit fitting formula [44], 
which is known to underestimate the nonlinear power on 
quasilinear scales by up to 4,2% at z = 0, 1 [43]. There- 
fore, we consider these results to be more than reasonable 
agreement, especially since there is no free parameter be- 
ing fit in making these predictions (see paper I). 

To sec the Vainshtcin effect in more detail it is useful 
to check the suppression of the brane extrinsic curvature 
represented by V^C. Note that if in Eq. (4) the nonlinear 
terms and non-local terms are neglected, we have the 
linearized relation 

Thus it is convenient to define 

_ arj 



(36) 



3(77-1) 

and study the suppression of c with respect to 5 expected 
from the Vainshtcin mechanism at small scales. In Fig- 
ure 7, we plot the ratio of the power spectrum of c, Pc, 
and the nonlinear density power spectrum P,n. Again we 
see that the ratio is 1 in small k limit, whereas it ap- 
proaches zero as k increases. The suppression is more 
appreciable at 2;=0 than z~\, as expected. 

Given the complexities of simulating the nonlinear 
PDE for C, Eq. (4), one may want to resort to approx- 
imations, and the spherical approximation immediately 
comes to mind. In this case one can convert Eq. (4) into 
a local equation in terms of V'^C at small scales. Indeed, 
spherical symmetry means that (VyC)^ = (1/3)(V^C)^, 
i.e. this corresponds to spherically average the Fourier 
space kernel in Eq. (29) over /i, giving 



2?^ — 1 



(l-/3V-i)<5, 



(37) 



which can be solved very easily at small scales (where 
nonlocal terms can be neglected), leading to the following 
modified Poisson equation [6, 8] 



3 rj 
2~ 



1 /3a 
2^VT 



(38) 




512 Mpc/h,z = l 
512 Mpc/h.z=0 
128 Mpc/h.z=l 
128 Mpc/;!,z=0 



k {h Mpc 

FIG. 7: The ratio of power spectrum of the normalized brane- 
bending mode c (see Eq. 36), Pc, to the nonhnear density 
power spectrum Pm. We show Pc from the 512 and 128 Mpc/h 
simulations at 2: = and 1. At small k the ratio is close to 1, 
but c is increasingly suppressed at high k. The suppression is 
also more pronounced for lower redshift, as expected. 
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FIG. 8: Ratio of PniDGP in the spherical approximation to 
the full solution as a function of scale for a simulation in a 
box of size 512/t~^Mpc (blue circles) and 128/i~^Mpc (red 
triangles) at z = 0. 



In practice we do include the nonlocal linear terms (al- 
though they are negligible for the box size we will show 
results for); then solving Eq. (38) is equivalent to sim- 
ply using the first iteration of our FFT-relaxation algo- 
rithm. Indeed, for our choice of w = 1/3, the first step 
corresponds to neglecting the quadrupole in the nonlin- 
ear nonlocal terms and solve for the monopole. This is 
similar to the method used in [45] to approximate the 
nonlinear Poisson equation for the DGP model and more 
generic degravitation models. 

Figure 8 shows the results of running our code in 
spherical mode compared to the full solution. Overall 
we sec good agreement, at approximately the 10% level. 
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However, there are important differences to note. First, 
there is disagreement between the two box sizes used 
(128/i~^Mpc and 512/i~^Mpc); second, the large-scale 
power of the two simulations differ among themselves and 
from the full solution. The reason for these results can 
be understood as follows. Making the spherical approx- 
imation corresponds to replacing (1 — /x^) —> 2/3 in the 
Fourier space kernel in Eq. (29), and that means there 
is considerable back reaction on linear modes from small 
scales. Indeed, defining 



$(fc) 



Ai']fc?C(ki)fc|C(k2), (39) 



where [S-q] = SuO^ — ^i — ^2), we are interested in esti- 
mating how much such a term would contribute to the 
growth factor at large scales. To do so, it is convenient 
to cross correlate $ with the large scale density pertur- 
bation, 



;$(fe),5(fe')) = Su{k + k') 



3(77-1) 



(40) 



where we used Eq. (35) to relate V^C to 6 at large scales, 
fx is now the cosine of the angle between k — q and q, and 
B denotes the bispectrum of density perturbations. Now 
as A; ^ 0, we have 



(1-M^ 



fc2 



(41) 



where is the cosine of the angle between k and q. 
Note that this geometric factor suppresses the contribu- 
tion of the nonlinear term to large scale modes. We are 
interested in the terms proportional to the large scale 
power spectrum, as we are looking effectively for the con- 
tribution in $(A;) proportional to 5{k). Using the second- 
order perturbative expression for the bispectrum we have 
that as fc — > such contribution reads 



($(5) cx /cV,2p(fc), 



(42) 



where ~ j qP [q) / <f . We thus see that in the full 
theory the contribution of these terms to large scales is 
suppressed by the factor fc^a^ at large scales, ensuring 
the validity of linear perturbation theory at sufficiently 
large scales. 

However, things change drastically if we take the spher- 
ical approximation (1 — ji^) — > 2/3, because the suppres- 
sion factor from Eq. (41) is no longer present, leading to 
instead 



( ^spherical (5 ) OC (T^F(fc), 



(43) 



where = J (fiqP{q). We thus have a contribution to 
large scale modes which is formally divergent, propor- 
tional to the variance of the density field and not sup- 
pressed by k. That's why the growth factor at large 
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FIG. 9; Dimensionless power spectrum of brane-bending met- 
ric perturbations Ave (see Eq. 44) as a function of scale at 
z — 0,1. Despite the fact that the dynamics of C is nonlinear, 
metric perturbations remain in the weak field limit. 



scales is more affected for the smaller simulation box 
(128 /i^^ Mpc) where there are smaller scale modes con- 
tributing to a larger cr^ than for the 512/i~^Mpc box. 
We also checked that the correction to the linear growth 
factor is smaller at higher redshift where is smaller 
than at 2 = 0. 

We conclude that the spherical approximation is rea- 
sonable at the 10% level but one should be careful that 
the solution is affected at all scales (including linear). 
From Fig. 8 we see that the Vainshtein effect is a bit 
stronger in the spherical approximation; this can be un- 
derstood from a calculation of the parameter Qcs in pa- 
per I, replacing the kernel (1 — /i^) by a constant increases 
Qoff ■ That's why despite having larger power than linear 
the spherical model approximates the full solution at the 
smallest scales (see Fig. 8). 

As the onset of the Vainshtein mechanism is related to 
nonlinearities in the dynamics of C it is worth checking to 
what extent the metric perturbations related to C arc in 
the weak field limit or not. As discussed in paper I, while 
nonlinearities in the dynamics of C are strong (mean- 
ing the extrinsic curvature fluctuations SK/K become 
larger than unity), the metric perturbations should re- 
main small throughout evolution. To check this, we con- 
struct the dimensionless power spectrum of the metric co- 
efficient (brane-bending mode) related to C, 54^ = VjC, 



(44) 



Figure 9 shows Aye as a function of k for z = 0, 1, con- 
firming that in fact Ave decreases as k increases. There 
are two reasons for this: first, the Vainshtein mechanism 
suppresses C in the nonlinear regime; second, extrinsic 
curvature fluctuations SK/K ^ {k/aH){k/a)C go non- 
linear at scales where (k/aH) ^ 1. From these numerical 
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ki = 0.1/iMpc~ and k2 = 0.2/iMpc~ as a function 
of angle 9 between fei and ^2. What we compare is the 
ratio of the reduced bispectra, defined as 



Q 



Bi 



23 



123 



iPlP2+P2P3 + P3Pl)' 



(45) 



where B denotes the bispectruni and subindices label 
wavectors, e.g. Pi = P{ki). The solid lines in Fig. 10 
show the expected ratio from perturbation theory as cal- 
culated in paper I (note however that fi^ = 0.27 here 
instead of 0.2), while symbols with error bars denote mea- 
surements in our A'^-body simulations (4 realizations with 
box size 1280 Mpc) for linear and nonlinear DGP 
models. The very good agreement is another (nontriv- 
ial) check of our numerical code. 



VI. DARK MATTER HALOS 



Mass Function 



FIG. 10: Ratio of bispectrum (nonlinear divided by lin- 
ear DGP) for triangles with k\ = 0.1/tMpc^^ and k2 — 
0.2 h Mpc^^ as a function of angle 6 between fci and fe2. Solid 
lines denote the perturbative prediction from paper I, symbols 
with error bars denote measurements in N-Body simulations. 



results we conclude that the Vainshtcin mechanism works 
as expected, and a cosmological background of perturba- 
tions is enough to suppress brane-bending mode fluctua- 
tions significantly at Mpc scales. 



V. THE BISPECTRUM 

The nonlinearities in the modified gravity sector (dy- 
namics of C) leave signatures in the observed non- 
Gaussian properties of density perturbations, as dis- 
cussed in detail in paper I. The prediction [8] is that the 
density field bispectrum should show an enhancement for 
isosceles triangles and no difference for squeezed triangles 
(where the nonlinear kernel in the modified gravity sec- 
tor, Eq. (29), vanishes since /i = ±1). 

Comparing the bispectrum in nonlinear DGP simula- 
tions versus standard gravity simulations is difficult be- 
cause they have different power spectrum amplitudes at 
low redshift, due to the difference in growth factors. As a 
result, the standard (i.e. present in GR) nonlinear effects 
will be slightly different and can mask the nonlinear ef- 
fects we are interested in here if those differences are not 
taken into account accurately enough. The best way to 
see the bispectrum enhancement is to compare nonlin- 
ear and linearized DGP simulations, which only differ in 
the presence of the nonlinearities in the modified gravity 
sector. 

Figure 10 shows such a comparison for triangles with 



We identify dark matter halos from the simulations 
using the Friend-Of-Fricnd algorithm [46] with linking 
length equal to 0.2 times the interparticle separation. 
We construct the mass function (dn/dlnM), the num- 
ber density of halos per unit comoving volume per unit 
interval of InM, for halos that have more than 20 parti- 
cles. We use box sizes of 128 ft."^ Mpc, 256 h~'^ Mpc and 
512 Mpc. We run 4 simulations for each box size and 
model. 

To better see the differences among the models we run, 
we divide our simulation results by the Sheth-Tormen 
(hereafter ST, [47]) mass function for the GR case 



dri: 



ST 



Pn 



dlnz 



dlnM MdlnM 



(46) 



where pm is the background density of matter, and 
V = [Sc/ctY, with 6c ~ 1-68 the linear overdensity above 
which a region will collapse according to the spherical 
collapse model, and cr^ is the variance of density fluctua- 
tions smoothed with a top-hat filter (with radius related 
to halo mass by M = (47r/3)pni-R'^)- The function vf{v) 
is a generalization of the standard Press-Schechter func- 
tion [49] motivated by the ellipsoidal collapse model and 
free parameters chosen to fit the numerical simulations, 
and is given by 



-v' 12 



(47) 



where v' = av, with a = 0.707. The numerical values of 
A and p are 0.322 and 0.3 respectively. 

The left panel of Fig. 11 shows the measured mass 
functions from the numerical simulations for the differ- 
ent models normalized as discussed above. To minimize 
the effects of poor mass resolution we use results from 
each box size for halo that contain more than about 300 
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FIG. 11: Left panel: mass functions measured in numerical simulations for the different models, as labeled, divided by a 
reference ST mass function. Right panel: the measured mass function for the three non standard models divided by the GR 
mass function from the simulations. The solid lines denote the predictions for the ST mass function ratios based on the 
renormalized linear spectra calculated in Paper I. 



particles; when more than one box contributes to a given 
bin in mass, we average results from different box sizes 
weighted by the volume of the corresponding boxes. Still, 
in the low mass regime (M < lO^^'^/i'^M©), the differ- 
ences for the GR case from the ST mass function arc most 
likely due to poor mass resolution, but again we expect 
the relative differences among different models (which 
are very small) to be an accurate representation of what 
would happen in higher resolution runs. At high mass 
(M > lO^^-^/i-^Afo), the enhancement for the GR case 
over the ST mass function is consistent with what is seen 
using Gadget (see e.g. [40]), and it is in this regime where 
the deviations among the different models is largest. 

The behavior seen at the high mass end is as expected, 
since the GR model has the largest as, followed by the 
GRH model, and then IDGP and nlDGP models that 
share the same linear normalization. The difference be- 
tween IDGP and nlDGP models is reasonable given the 
Vainshtein effect, which makes the power spectrum of the 
nlDGP case larger than IDGP. At low mass, the differ- 
ences are much smaller, similar to what happens with 
GR models of different cosmological parameters where 
the mass function does not evolve much for M <C M*. 
In this regime the nlDGP model has a larger mass func- 
tion than the IDGP model, again this is expected from 
the Vainshtein mechanism, and indeed the nlDGP model 
agrees with the GRH model at small mass, as expected 
from the matching of their power spectra at small scales. 

Motivated by these results, the right panel of Fig. 11 



shows the mass function calculated from Eq. (46) using 
the renormalized linear power spectra as calculated in 
paper I. Since 6c = 1.68 works very well for different cos- 
mologies in the case of GR (e.g. [48]), we take this value 
for all models we consider regardless whether gravity is 
modified or not. Corrections to this in the DGP model 
are expected to be small [6]. While we do not expect 
this mapping to work in too much detail, the idea is that 
we are including the Vainshtein mechanism in the renor- 
malized linear spectra for nlDGP derived in paper I and 
using the ellipsoidal collapse model assumed in the ST 
mass function to take into account the standard grav- 
ity nonlinearities. For work on spherical collapse model 
implications for the mass function in modified gravity 
models see [50]. The predictions of paper I do reflect 
roughly the features seen in the measurements. Indeed, 
the relative ordering of the mass function measurements 
agree with expectations based on the predictions. 



B. Halo Bias 

We now discuss the difference in the large-scale bias 
of dark matter halos for the different models discussed 
above. Due to our relatively low particle number our 
mass resolution is limited, therefore we simply concen- 
trate on all halos with more than 20 particles to com- 
pute the bias. The halo bias is obtained from the ratio 
of the cross halo-matter power spectrum and the dark 
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FIG. 12: The halo bias for halos with mass M > 1.5 x 
IO^^Mq/Zi (top four curves) and M > 1.8 x IO^^A/q/Zi (bot- 
tom four curves) for different models. 

matter power spectrum, to avoid accounting from pos- 
sibly nontrivial shot noise in the halo-halo power spec- 
trum. Figure 12 shows the halo bias from simulations of 
box size 256/i~^Mpc and 128 Mpc, with the same 
threshold of 20 particles per halo this corresponds to 
b{M > 1.5 X lO^^Me/h) and b{M > 1.8 x lO^Mg/Zi), 
respectively. The error bars correspond to the mean over 
four realizations. 

At fixed mass threshold, the relative magnitude of 
the bias agrees with expectations based on the peak- 
background split [51, 52], which says that the bias is 
proportional to the slope of the mass function. One can 
check that the ordering of the bias factors in Fig. 12 at 
fixed threshold mass qualitatively agrees with this simple 
prescription by looking at the slopes of the mass functions 
in the left and right panels of Fig. 11. 

In [53] it is pointed out that a scale dependence of the 
growth factor (as it happens in DGP models at large 
scales due to nonlocal terms) can give rise to a scale 
dependence in the halo (or galaxy) bias. Assuming no 
merging and that halos move with unbiased velocities, 
the halo bias b{z,'k) in Fourier space evolves as [53] 

where bo denotes the halo bias at formation time, and 
Sni is the density contrast of the matter. Even if bo is 
independent of scale, the ratio of matter density fluctu- 
ations at different times is expected to show some scale 
dependence due to nonlocal terms in the modified Pois- 
son equation (see paper I). However, we were unable to 



test this prediction because 1) the scale dependence in 
the growth factor is small and kicks in at large scales, 
2) to probe large scales, large boxes are required, and 
our limits on computational resources required means we 
cannot use large number of particles, which leads to poor 
mass resolution, 3) for large halo masses (poor mass res- 
olution), halos form late and there is little time between 
formation and z = to see a relative scale dependence 
induced by the ratio of matter fluctuations in Eq. (48). 
As a result of these factors, the predicted scale depen- 
dence is too small and consistent with the uncertainties 
from our simulations. 



VII. CONCLUSIONS 

We developed a numerical algorithm, which we call 
FFT-relaxation, to solve the nonlinear field equations 
(Eqs. 3-4) in the DGP model. This enabled us to run fully 
nonlinear iV-body simulations of this model and compare 
to simulations where the field equations are linearized, 
where standard gravity is used only for perturbations 
(but not for the background), and to fully standard grav- 
ity simulations. Our numerical algorithm includes not 
only the new nonlinear couplings in the modified grav- 
ity sector but also linear but nonlocal terms that lead to 
5D behavior at large scales in the quasistatic approxima- 
tion. We expect that our implementation may be useful 
in other scenarios of brane-induced gravity or massive- 
gravity related models, e.g. [34-36], which have both of 
these features. We extracted the power spectrum, bis- 
pectrum, mass function and halo bias from the iV-body 
simulations. 

From the simulation results, we confirm the expecta- 
tions based on analytic arguments that the Vainshtein 
mechanism does operate as anticipated, with the den- 
sity power spectrum approaching that of standard gravity 
within a modified background evolution in the nonlinear 
regime. The transition is very broad and there is no well 
defined Vainshtein scale, but roughly this corresponds 
to wavenumbers fc* ~ 2/iMpc~^ at redshift z ~ 1 and 
fc* ~ 1 /i Mpc~^ at z = 0. We checked that while extrinsic 
curvature fluctuations go nonlinear, and the dynamics of 
the brane-bending mode C receives important nonlinear 
corrections, this extra scalar mode does get suppressed 
compared to density perturbations, effectively decoupling 
from the standard gravity sector. At the same time, there 
is no violation of the weak fleld limit for metric perturba- 
tions associated with C. We also quantified the use of the 
spherical approximation to the nonlinear dynamics of C 
(see e.g. [45] for application of this to more generic mod- 
els), and showed that the results on the density power 
spectrum show differences of order 10% at all scales, de- 
pending on box size and resolution. The differences are 
not restricted to small scales as the approximation intro- 
duces a spurious enhanced coupling of large-scale modes 
to small scale modes. 

In the nonlinear regime, modified gravity models be- 
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come difficult to distinguish from standard gravity for 
a fixed expansion history, and therefore this is not the 
best place to look for modifications of gravity in obser- 
vations. The best way to test these models is on weakly 
nonlinear scales, we showed that the bispectrum in our 
simulation presents a small signature of modified gravity 
due to the additional nonlinear couplings, as predicted 
by perturbation theory in paper I; however, this is too 
small to be tested in present observations. We also stud- 
ied the abundance and clustering of dark matter halos 
and found that the largest differences in the mass func- 
tion are seen at large masses, as anticipated from the 
results in the power spectrum. The differences seen in 
the halo bias are consistent with expectations based on 
the peak-background split. 

We tested the calculations presented in paper I that 
rely on a renormalization of the linear spectrum due to 
nonlinear effects in the modified gravity sector, which 
are induced by the Vainshtein mechanism and modify 
the effective gravitational constant from linear theory. 
The predictions for the nonlinear spectrum are in very 
reasonable agreement with our simulations. A similar 
prediction for the mass function shows the right trends 
although our tests are modest due to limited simulation 
volume and mass resolution. We hope to improve on 
these tests to sharpen theoretical predictions of modified 
gravity in the near future. 

Our results suggest that the framework laid down in 
paper I, where the nonlinear effects of gravity beyond 
general relativity arc included in the running of the grav- 
itational constant, is quite accurate. As discussed in pa- 
per I, this is in fact a result of the fact that modifica- 
tions beyond such running (such as extra contributions 
to non-Gaussianity) are higher-order in SG/G and are 
highly suppressed at small scales. This general argument 
should be applicable to other theories of gravity, leading 
to a powerful way to test gravity theories against obser- 
vations. 



While this paper and its companion paper I [8] were 
slowly written down, preprint [54] appeared on arXiv, 
where N-body simulations of the DGP model arc also 
presented. His results on the power spectrum and halo 
mass function are broadly consistent with ours. 
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APPENDIX A: OTHER ASPECTS OF THE 
FIELD EQUATIONS 

In this appendix, we would like to compare the FFT- 
relaxation method with more common relaxation meth- 
ods. Since the difhculty in solving the Eq. (4) stems from 
the nonlinear terms and it is difficult to implement the 
nonlocal terms in the usual finite differencing schemes, in 
this Appendix we discard the nonlocal terms for simplic- 
ity. Thus here we are interested in solving the simplified 
equation 

(V^C)^ + aV^G - (V,,C)2 = fcilj. (Al) 



1. Ellipticity of the Equation 

Before attempting to solve Eq. (Al), some comments 
are in order. First, Eq. (Al) belongs to the class of fully 
nonlinear partial differential equations (PDF), since it is 
nonlinear in the highest derivative of the unknown. Fully 
nonlinear PDEs are usually more difficult to deal with 
than those nonlinear in the field, such as the field equa- 
tion in f{R) models [13] or nonlinear in lower derivatives 
of the field. So far there is a general solution theory for 
this kind of fully nonlinear PDEs only if Eq. (Al) is el- 
liptic in the whole domain [55]; however, to determine 
whether Eq. (Al) is elliptic or not is nontrivial. 

In general for a PDF F{G) = 0, in our case 

F{G) = (V^C)^ + aV^C - (V,,C)2 - ^^^—^1,5, (A2) 

one needs to solve for the eigenvalues of the matrix 
^v,,c Pv^yC ^v,,c\ 

PVy^C PVyyC y , C ] , (^3) 

^v,.,c Pv,^yC Pv,,c) 

where i^y q denotes partial differentiation of F with 
respect to V xxG , etc. The equation F{G) = is elliptic 
if all the eigenvalues are of the same sign, hyperbolic if 
the sign of one of the eigenvalue differs from the rest, and 
parabolic if one of the eigenvalues is zero. 

However, in 3D, the resultant equation depends on 
the field C, and we cannot determine its nature of the 
equation without knowing C! Nonetheless, we can check 
the ellipticity of the equation using the solution we ob- 
tained from the FFT-relaxation method. We do this 
test in a 64 Mpc simulation box using 256^ particles. 
We input the convergent solution C(x) into the matrix 
Eq. (A3) and solve for its eigenvalues at every grid point. 
We then check if the signs of the eigenvalues arc the same 



16 



or not. Wc find that the equation is elhptic at every point 
in the simulation box. 

It is important to note that violation of cUipticity hap- 
pens when the bad fraction is non-zero. In our test we 
used = 1/3 and therefore the bad faction is zero all 
the time. If = is used instead, the bad fraction will 
be non-zero as shown in the upper panel of Fig. 4, and 
violation of ellipticity (by substantial fraction of points, 
of order 10%) will occur in the same time interval during 
which the bad fraction is non-zero. Thus satisfying the 
ellipticity condition is another self-consistent check of the 
FFT-relaxation method. 

If Eq. (Al) were in 2D rather than 3D, its nonlinear 
structure would be the same as the well-known Monge- 
Ampere equation [56]. This simplification is possible if 
we artificially assume that 6 is constant along the z direc- 
tion, with C a function of x and y alone. In 2D. Eq. (Al) 
reduces to 

2{V^^CVyyC-iV^yC)^)+a{V^^C + VyyC) = ^^^^^ S 

and the determinant of the corresponding matrix 
Eq. (A3) in 2D, after making use of Eq. (A4), is 

a' + fcll<S. (A5) 
V 

One can check that Eq. (A5) is always positive even for 
S — —1. This shows that the eigenvalues are of the same 
sign, and thus Eq. (A4) is elliptic. 

2. On Gauss-Seidel Relaxation 

Now we turn to solving Eq. (Al) using a finite differ- 
ence method. One of the potential methods for solving 



nonlinear PDEs is the Newton-Gauss-Seidel-relaxation 
method [57], which for our case it reads 



where C"j^ denotes the value of C at the grid point 
{i,j,k) in the n*^ iteration and F'{C^^) represents the 
derivative of F with respect to C^jj.- We find that this 
method diverges very quickly. This is probably due to 
the fact that upon discretization, Eq. (Al) is a quadratic 
equation in the discrete variable Cijk, which has two, one 
or zero real roots. If the starting trial solution is not good 
enough, it is impossible for different points to converge 
to the same branch of solution. 

Recently, algorithms were developed by mathemati- 
cians [58] to attack the Monge-Ampere equation, but 
most of them are difficult to implement, albeit more rig- 
orous than the method we developed in this paper. How- 
ever, equations similar to Eq. (A4) also arise in other 
fields such as meteorology, where a 2D equation called 
the "balance equation" [59, 60] used to model the wind 
flow is of Monge-Ampere type. The approach followed 
in those papers is simple enough to implement. They 
used Gauss-Seidel relaxation together with the quadratic 
formula instead of Newton's method. This chooses the 
branch of the solution explicitly and has better control 
over which branch the solution converges to. In what fol- 
lows, we will discretize Eq. (Al) similarly to how is done 
in [59] and generalize it to 3D. 

Le us now discuss the details of the calculations. We 
first write. 



(V2C)2-(VyC)' 

= 2[V^^CVyyC + V,,CV,,C + VyyCV,,C] " 2[{V^yCf + {V + [V y ,Cf] 

= (V,,C - \JyyCf + (V,,C - V,,C)2 +a^- {\JyyC - V,,C)2] - 2[(V,yC)2 + (V,,C)2 + {\Jy,C)\ 

where in the last line we have denoted /i = xxC -f Vj^j,C, u = xxC + V^^C, and a ~ ^ yyC + V zzC respectively. 
Now Eq. (A2) can be expressed as 

a 1 1 - - 

F = -(a* + + a) + -(m' + + - -[(V,,C - ^yyCf + (V,,C - ^zzCf + (V^^C - \/ z^Cf] 



i{r,-l) 



V 



, a. 9 , a. 9 , a. 9 3a 
+ 2 )' + (^^ + 2 ) + ('^ + 2 ) — 



- 2 - VyyCf + (V,,C - VzzCf + {VyyC ~ V.^Cf] 



~2[{VxyCf + (V,,C)2 + [VyzCf] 



V 



(A7) 



The advantage of writing F in this way is that upon discretization, the central term Cijk only appears in the first 
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square bracket. Suppose we write the difference equation as a quadratic equation in Cj"j,, and solve for its root to 
get the new C^""^^ (note that we assume that the iteration is Jacobi-hke, i.e. wc use only the old C" to update the 
central term C^jf^; we will relax this assumption later). Let's denote the residual before the n**^ iteration be F", and 
the residual after the n"^ iteration be Explicitly, we have 



2 



F 



n+l 



(A8) 
(A9) 



where we have denoted the non-central terms Z{C'^). which arc unaltered in the iteration. Subtracting F^ from 
we have 



(AlO) 



Next, we discretize the variables and the derivative operators, and write Cj"^^ as CJJj. + ^^"j.. Then, for example 



■'ijk 



we can express 



n+1 _ ^t+l.jfc 
f^ijk - 



a 



a 



i,j + l,k 



/l2 



(All) 



where h = ai//i and h is the size of the grid. The extra factor aH originates from our definition of V as V/{aH). 
Thus we arrive at a quadratic equation in f"-. 



6 



3a 



ijk 



'ijk 



'ijk 



ijk 



/)4 

24 ^ 



Fr,k)=0. 



(A12) 



We cannot solve this equation since F[^^^ is unknown. The next approximation we make is to set -f/J^^ to zero. 
This is reasonable as we assume the next iteration will give a smaller residual than the previous one. Thus the final 
equation we get is 



i^Tik)"^ — TiP'i'ik + Kik + Kjk + -Tr)^iik + TTjPijk — 



(A13) 



From Eq. (AIS), we can solve for ^^j. and then get new C"jl^ . However, the quadratic nature of the equation 
introduces another complication — which root to take. Obviously, when we take the limit i^^'j^. —^ 0, we should take 
the root such that ^^j. — *■ 0. Thus we shall take the plus sign of the square root. Again, we may be plagued by 
bad points since the discriminant may not always be non-negative. Rearranging the terms as in Sec. Ill C may help 
resolve this issue. We will not, however, pursue further investigation since we will see that this method docs not quite 
converge even when the bad fraction is zero. 



So far, we used a Jacobi approach, i.e. in updating 
C"+^ we only use old values of C". Now let's consider 
a Gauss-Seidel approach, in which we will use the new 
(7"+^ whenever they are available during the iteration. 
We find that this can improve the convergence rate sig- 
nificantly, as it is normally the case. 

We compare the convergence of the GS-relaxation and 
the FFT-relaxation method in Fig. 13, where we plot 
the normalized residual versus time taken for these two 
methods. We start from the same initial linear solution. 
We perform this test on a box of 64 Mpc at z = with 
w = 1/3 for FFT-relaxation. In both methods, there are 
no bad points (in fact for the FFT-relaxation method. 



a small bad fraction appears in the first few cycles, and 
then disappears). 

We see that although GS-relaxation takes less time to 
run in each cycle, its normalized residual settles to some 
value of order 1. In fact the normalized residual rises very 
slowly as we continue to iterate. On the other hand, for 
the FFT-relaxation method (although takes longer time 
for each cycle) the residual keeps on decreasing, especially 
after the first few cycles. During our simulations, except 
for the very first time step, we use the previous time 
step solution as the trial solution for the FFT-relaxation 
method, which is much better than the linear solution 
since we then are in the steep-slope regime in Fig. 13. 
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— GS-relax 
--■ FFT-relax 




10^ 

time (sec) 

FIG. 13: Normalized residuals as a function of time for the 
FFT-relaxation and GS-relaxation methods. For the GS- 
relaxation method the normalized residual saturates at some 
value of order 1, in this case about 1.4, while for the FFT- 
relaxation method the residual keeps on decreasing. Note that 
in practice, we work on the steep part of the FFT-relax curve 
when we use the previous time step solution to initialize the 
relaxation method. 



However, it may be useful to combine both methods to 
achieve greater convergence speed. 



It's not clear to us why the GS-relaxation method 
together with the quadratic formula seems to work 
in [59, 60], although there are significant differences with 
our application here. First, the balance equation is 2D 
rather than 3D. Secondly, the number of grid points con- 
sidered in [59, 60] are only about 10'^, so the resolution is 
rather low. Last, and perhaps most importantly, the per- 
formance of these nonlinear relaxation methods depend 
rather sensitively on the values of the parameters such 
as a and 3(?7 — 1)?7~^ and how rough the data S is. In 
cosmology, 6 is not very regular. To summarize, we pre- 
sented a simple discretization scheme for GS-relaxation 
method, and it fails to solve Eq. (Al) satisfactorily. How- 
ever, given the complexity of Eq. (Al), we don't rule 
out the possibility that a more thoughtful discretization 
scheme would work. In fact, while this paper was be- 
ing finished, preprint [54] appeared on the arXiv showing 
that a convergent method based on Newton's method can 
be constructed for these field equations based on multi- 
grid techniques. 
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